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Abstract. Within investigating the multiple charge spreading generalizing the 
Bertaut approach, a set of confined spreading functions with a polynomial 
behaviour, but defined so as to enhance the rate of convergence of Coulomb series 
even upon a single spreading, is proposed. It is shown that multiple spreading 
is ultimately effective especially in the case when the spreading functions of 
neighbouring point charges overlap. In the cases of a simple exponential and 
a Gaussian spreading functions the effect of multiplicity of spreading on the rate 
of convergence is discussed along with an additional optimization of the spreading 
parameter in dependence on the cut-off parameters of lattice summation. All the 
effects are demonstrated on a simple model NaCl structure. 



PACS numbers: 02.30.Lt, 02.30.Uu, 61.50.Ah, 61.50.Lt 
1. Introduction 

It is known that the Ewald approach pQ is the most widespread implementation of 
the Poisson summation formula in crystals [21 [3l [21 [5] . Bearing in mind that such 
a treatment proposes that the overall summation is divided into sums over direct 
and reciprocal space, the enhancement of the rate of convergence is traditional in this 
problem and assumes the investigation of ranges of summation in either of those sums. 
Indeed, Epstein [B] discussing this problem for the first time has proposed that both 
the sums must be separated by a dimensionless parameter unity, in agreement with 
conventional mathematical approaches. According to Ewald pQ, the corresponding 
parameter of splitting can still be chosen as a unique, but variable, so as to provide the 
most rapid rate of convergence of both the sums. This proposal remained appreciable 
for a long time [7|. However, in the last years one more standpoint arises, keeping 
in mind that the splitting parameter can also depend on the cut-off parameters of 
both the summations in question [SI IS]- I n particular, such a treatment can reduce 
the computational efforts associated with the dimension of a supercell employed in 
molecular dynamics [10l E] • 

It is important that the Ewald approach can be regarded as an effect of charge 
spreading with a Gaussian spreading function [2 El [11]. In this connection, the Bertaut 
treatment |12j is obviously the extension of the Ewald one to an arbitrary spreading 
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function. It implies that the problems of adjustable parameters of splitting still remain 
actual in this generalized treatment as well [HI H3] • Of course, these problems become 
immaterial if the spreading functions applied to different point charges in a point- 
charge lattice do not overlap [TU Q3] . In this case the sum over direct space is absent 
and so the error arising upon truncating the summation over reciprocal space is the 
only subject of interest [T5 ] [TB I \T7 \ fT8 l fT9] . 

It is worth noting that the charge spreading with a certain spreading function 
applied to all couples of interacting charges in the lattice was originally discussed 
[T2l 03]. On the other hand, it turns out that spreading the charges generating 
the potential field is sufficient for enhancing the rate of convergence of the lattice 
series [T3j [20j [21]. Nevertheless, the application of spreading to all the charges in 
the expression for the Coulomb energy appears to be somewhat more efficient |14j . 
The reason of this efficiency can be understood from the fact that, by symmetry, the 
latter effect may be treated as a double spreading of charges generating the potential 
field. As a result, the idea of a multiple charge spreading arises as a next step towards 
achieving the faster convergence [22] . 

In the present paper the problem of enhancing the rate of convergence is 
discussed in detail. We consider different principal classes of spreading functions, 
with concentrating our attention on the effects of multiple spreading. As particular 
cases of spreading functions extended to infinity, here we discuss a Gaussian spreading 
function and a simple exponential spreading. Furthermore, we examine different types 
of confined spreading functions in a sequence of enhancing their convergence efficiency, 
providing that the possibility of their overlap is accessible. To our mind, the latter 
effect extends the original ideas of Bertaut Q2] [14] . 

2. Basic relations describing the multiple charge spreading 

For convenience, here we compile some results describing the multiple charge spreading 
and obtained earlier [22]. Every perfect crystal can be specified by a unit-cell charge 
distribution p(r) subject to the condition of neutrality of a unit cell : 



Jv 

where the integral is over the volume occupied by p(r). The corresponding structure 
factor as a function of a reciprocal lattice vector h is of the form 



where v is the unit cell volume. Then the relation F(h = 0) = readily follows from 
(TTJ) . In the particular case of a point-charge lattice the charge distribution is converted 
into 



where j runs over point charges qj belonging to a unit cell and located at positions 
bj, 5(r) is the Dirac delta function. 

The effect of charge spreading will be described by a spherically symmetric 
spreading function ad?"!) normalized by the condition 




(1) 




(2) 




(3) 




(4) 
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Here r — \r\. The application of <r(|r|) to the charge distribution in a unit cell results 
in an additional multiplier 

S{h) = J a(\r\) exp(~2TTihr) dr (5) 

modifying structure factor J5} as follows: 

F(h) -> F(h)S n (h) (6) 

if the spreading of interest is performed n times in a consecutive manner. 

The corresponding potential effect of a multiple spreading in direct space is 
described by the function 

11 u J |Jl + ri + ... + r n | V ; 

As a result, it is found that the electrostatic potential field arising within such a 
procedure of spreading can be written as 



rr / X 1 F(h)S n (k) , , . 



+ E7^n)^»-{^^(0)} b , (8) 

where the prime on the summation sigh over reciprocal lattice vectors h implies missing 
the term at ft. = 0, the parameter i runs over vectors Ri determining the Bravais lattice 
and specifying different unit cells, we introduce the notation Ri = Ri + r\ — r, the 
prime on the summation sigh over i means that a possible singularity associated with 
the denominator in the summand must be omitted as well, 

= (R), (9) 

the last term on the right-hand side of ([8]) is the correction associated with a point 
charge qj if it happens at r. Taking potential ([5]) into account, we write down the 
Coulomb energy per unit cell in the form 

«wgh)fS^) lwf 

x ^-«Ee do) 

\Ril2\ Z j 

where the summation over j is over all point charges qj in the unit cell and Rn 2 = 
Ri + r x - r 2 . 

The recurrence relations associated with n^(R) and £?( n )(0), respectively, take 
the form 

cy poo rR+r 
Q{n) {R) = _ a{r)rdr Q^- 1 ){y)ydy, (11) 

M JO J\R-r\ 

/>Q<Q 

n {n) (0) = 4n tT(r)tJ<- n - 1 \r)r i dr. (12) 
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The case of n — 1 is then straightforward [22J and is described by 

/>oo 

W^(R)=4w / a(r)r(r - R)dr, 
Jr. 

/>oo 

n (1) (0) = 4vr / a{r)rdr, 
Jo 

providing that Q^(R) = l/R. In the case of n — 2 one can obtain 



(13) 
(14) 



W (2) (R) = 4tt 2 



dri / dr 2 A(r 1 ,r 2 ) 



dri 



R-n 



dr 2 



x A( ri ,r 2 ) 



dri 



dr 2 B(r 1 ,r 2 ) 



R+ri 



n (2) {0) =32vr 2 / a{r x )rxdrx \ a{r 2 ){r 2 ) 2 dr 



where the following definitions 

A(rx,r 2 ) = a(rx)a(r 2 )rir 2 (R - n - r 2 ) 2 , 
£(n,r 2 ) = cr(ri)cr(r2)rir2(i? + r x - r 2 ) 2 
are taken into account and notation ([9]) is also utilized. 

Employed to the simple exponential spreading specified by 
a 3 

CT ( r ) = ^ ex P(~ ar )' 
the latter results give rise to 

S(h) 







1 + 





wW(z)= (l + |)exp(-z) 



H ) exp(— z), 

16 16 48/ v ; 



(15) 
(16) 

(17) 
(18) 

(19) 

(20) 
(21) 

(22) 



T ~_/o\ . . / 193z 65z '6lz' J z~* z" \ . . 

W (3) z = H h 1 1 1 exp (-z), 23 

K ' V 256 256 768 192 3840/ PV 7 ' v ; 



37z 



fl«(0) = | 



fl«(0)-g. ^(o) = ||. 



(24) 



Here we go over to a new dimensionless variable 2 = aR, so that the notations 
introduced above are modified as follows: 

W {n) (R) = W (n) {z). (25) 
The other important case is associated with a Gaussian spreading function 

1 2 \3/ 2 

(26) 



M r ) = ( 

\ 7T 



expl — nfi 2 r 2 ) 



The results appropriate to this case are of the form 

ir 2 \h\ 2 > 



S%{h) = exp(- 
W [n \z) = erfc(z) 



exp(— u 2 ) du, 



/}( n >(0) 



2/i 



(27) 
(28) 
(29) 



The fact that the final results are independent of n is the peculiar feature of the Ewald 
approach, as stressed earlier. 
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3. Spatially confined spreading 

Let us now consider spherically symmetric spreading functions bounded by a radius 
Ro'- 

( s f 9k( s )(r) at r < R , 
<W = < (30) 
10 at r > Rq, 

where we restrict ourselves to polynomials gk( s ) °f order k which in turn result in 
Sk( S )(h) oc |/i|~ s - The meaning of the parameter s is defined therefrom. Then it is 
evident that relation © contributing to still takes form (|2"5)) . but at z = R/Ro. 
Furthermore, it will appear that Wj. n \z) = as z > n. This fact is natural for 
the Coulomb potential generated by a spherically symmetric charge distribution and 
implies that the sum over i in (|8|) is actually finite and includes only unit cells nearest 
to a reference point . 

There are different polynomials discussed in the literature [12j [13l HH [16j HB1 EQl 
|2"T1 12"31 |2"4") . Interested in principal particular cases of spreading (fBT))) , we begin with a 
uniform spreading normalized by condition ([4]), as proposed by Bertaut [12] : 

5 °( 2 > {r) = 4^r (31) 

According to (0, relations (|3T))) and (pTTj) yield 

S 0{ 2)(h)^^(^-cosY), (32) 

where Y = 2Tr\h\R and so s — 2 herein. If n = 1, then the other quantities of interest, 
which are specified by the subscript k(s) and by the superscript (n), are as follows: 

(l-z) 2 (2 + z) 



at < z < 1, 



at 2 > 1, 



<i(0) = ^. (34) 

The fact that So(h) contains the factor enhances the rate of convergence of a 
series over h in ([8]) even at n = 1. This effect becomes stronger as n increases. In the 
original approach of Bertaut [12] the case of n = 2 is considered. Making use of (fl~5|) 
and (fT6|) . one can show that 

f (2-z) 4 (10 + 8z + z 2 ) 

<»,(,)= ^160 1 « t0 ^ S2 ' (35) 

[ at z > 2, 

<»,(») = (36) 

in this case. The event of n = 3 appears to be more complicated, though it is still 
described by formulae (fTTj) and (fT2"]) based on relations and ([3"5|) . The results can 
be cast in the form 

flgg( ) = J^° (37) 
°( 2 ) v ; 1280i? 

Mi(z) + M 2 (z) at < z < 1 , 
^o^)^) = "I M i( z ) at 1 < z < 3, (38) 



at z > 3, 
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where 

(3 -zf(72 + 81, + 18* + z 3) 
n y 53760 ' 1 ; 

M 2 ( Z ) = (l^) 6 (424 + 87 Z -6^-z3) 

v ; 17920 1 j 

It is important that, apart from enlarging the value of n, another way of enhancing 
the rate of convergence of the series over reciprocal vectors can be arrived at by 
increasing the order k of g^u) so as to enhance the degree s of Y in the denominator 
°f Su s )(h) [21] • In particular, one can show that 

3(1 - x) . . 

9 m (r)= 3 (41) 

where x = r/i?o, leads to 

*»W-p( J£ T !2 -*i'). W 



In this case 



z) 3 (l + z) at < z < 1, 
at z > 1, 



— ' " u f:^ (43) 



<3)(0) = ^- (44) 



Extending the consideration to the case of n — 2 here, we derive 

? (2) 
'l(3) 



«(0) = Jr, (45) 




M 4 (z) at0<z<l 
at 1 < z 
at z > 2, 



= <( M 3 («) at 1 < z < 2, (46) 



where 



Ma(2) _ (2-»)°(2 + 4, + ^ (47) 

M 4W . P + (48) 

Interested in polynomials of the lowest degree, the next several polynomials for 
the spreading function within this set are shown in |Appcndix A| Note that the couples 
of relations (J3TJ) and ([32]). (|4T|) and flU} and finally (jA.lj) and (|A.2|I from | Appendix A| 
naturally agree with results proposed in [131 HH EH [211 123] • One can see that the set 
of spreading functions is reduced to terms of the form (1 — x) k [T3l [TH [20] only at 
k = 1 and k — 2. This fact accounts for the known statement [16] that the simple 
form (1 — x) k is not yet efficient at k > 2. 

Note that various polynomials are possible for a given s if they are not restricted 
to the lowest degree. Of different polynomials associated with s = 3 and so competing 
with ([4T]) , here we consider the only one that is expected to be very effective [24] and 
is of the form 
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Its Fourier transform is also known: 

c /l\ 20/6sinT 2 + 4COSF . n 

Sa(3)W = ys(^2 y S1IlF J- ( 5 °) 

All the other quantities appropriate to the case can be obtained in the manner 
considered above. The case of n = 1 is described by 

f (l-z) 3 (3 + 4z + 3z 2 ) 

^1 1 a,0S = S1 ' (51) 

[ at 2 > 1, 

In the case of n = 2, we reach 

' M b (z) - M 6 (z) at0<2<l, 
W^l } (z) = \ M 5 (z) at 1 < 2 < 2, (54) 



where 



at 2 > 2, 



M 5 (,)^ 2 -^ 28 + 24z + 42z2 + lfe3 + 3 ^, (55) 
y ' 1512 ' v ; 

M.W . 5 "- j)T ^ 4 ' + j2 ' . (56) 



4. Optimization of spreading parameters 

It is important that every spreading function is characterized by a certain parameter. 
The problem of its optimization is traditional. In particular, the value of fi = 2^/ir/d 
describing a Gaussian spreading function was used as a unique for the NaCl structure, 
where d is the lattice spacing [UE]- The situation associated with a simple exponential 
spreading appears to be quite similar [T3J [2S1 US] • 

First of all, here we develop the treatment appropriate to this case. As mentioned 
earlier 5 , it is based on the Coulomb characteristic of a Bravais lattice, the parameter 
put forward by Harris and Monkhorst [57]. In order to discuss this approach, 
we consider a Bravais lattice composed of unit point charges and immersed in a 
neutralizing uniform background. It is easy to show that the interaction of a 
background with the bulk potential field vanishes and the same is right for the 
background contribution to the sum over reciprocal lattice in expression (p~0|) [5, 22, 27J. 
As a result, the effect of lattice summation in (jTOJ) is associated with the contribution 
of point charges alone. However, there is a remainder constituted of two simple finite 
terms there. One of those terms is determined by the last optional term in l|10[) . The 
other one is the contribution of a background to the integral over real space and so is 
described by the negative of the quantity 

o 2 poo 

GW = — / r 4 a(r)dr, (57) 
3w Jo 

keeping in mind that a background is of negative charge |22j . 

In the case at hand equation IJ2J) leads to F(h) = 1. It implies that both the series 
over reciprocal and direct lattices in the original energy expressions are divergent and 
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the spreading length, denoted in general as A, actually forms cut-off parameters for 
both of those series in formula (fTU)) . Therefore it is expedient to cast the expression 
for the Coulomb characteristic C of a Bravais lattice in the following schematic form: 



A „ 
— + BX 
A 



d n) +/2™(0) + C, 



(«), 



(58) 



where the terms on the left-hand side stand for the contributions of series over 
reciprocal and direct lattices, respectively, with singling out their principal dependence 
on A. Note that C is the value of the modulated lattice potential at the site of a unit 
point charge and so C/2 = £, where £ is the total energy per point charge in question 
[5]. It is significant that except for C, all the terms in relation ([55]) are positive. 
Moreover, according to definitions , (fl~2|) and (f57|) . one can readily show that 



Therefore G^' l) and 0^(0) may be treated as counterparts to the direct-lattice series 



Gi n) «A 2 , 
,(«) 



4 n) (0) oc A" 



(59) 



and to the reciprocal-lattice one in (|58p . respectively. Bearing in mind that C is a 
constant, we draw a conclusion that an optimum relation between the terms on the 
left-hand side of formula (|58[) is associated with a minimum of the right-hand side 
that is determined by the condition 

"' """ 12™ (0) =0. (60) 



dX 



G 



(n) 



The Gaussian spreading 

7T 



G u 



substituted into formula (|5T|) at a given n results in 

(61) 



Note that this functional form is still independent of n. In this case X — 1/ fi. Inserting 
relations and (pTjl into ([SO]) , we derive the result discussed earlier [H [5] : 

m °p' = rfh>-V 3 , ( 62 ) 

providing that the differentiation immediately with respect to /i in ([60]) is also 
admissible. 

In the cases of simple exponential spreading discussed above, the substitution of 
l[T9|) into (52|) yields 

87m 
^2 ' 



(63) 



Now we substitute issue (163)) along with one of the relations described by (|24j) into 
(|60p . The peculiar feature of relations l|55|) implies that both the differentiation with 
respect to a and the differentiation with respect to A = 1/a in (p0|) also lead to the 
same result: 



4tt\1/3 



2 (t) 
4- 



1/3 



16 



7T 

21t; 



1/3 



at n — 1 



at n 



at n = 3 



(64) 



Note that the value of a opt increases when n changes from 1 to 3. It means that 
the initial spreading function becomes more compact and thus the enhancing effect of 
spreading at n — 2 and further at n = 3 turns out to be somewhat restricted. 
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Another case of optimization arises if the values of the spreading parameter a, 
for definiteness, are supposed to depend on the cut-off parameter of summation m. If 
an integer value of m is common to both summations over reciprocal and direct space, 
then the result for the energy can be written as . Apart from this quantity, it is 
also advantageous to consider the value of £|™+ ma ] ^ where m a is a certain integer as 
well. The value of associated with a given m can be determined by the condition 

£ [rn+ ma] _ g[m] = Q (65) 

In order to understand this relation, one should point out that at least two different 
situations arise here. Indeed, if m a = 1, then formula (165j) is apparently reduced to 
the condition of stability with respect to small variations of m. The opposite limiting 
case arises when ui a ^> 1. This case corresponds to the fact that ((65)) evaluates all 
the rest removed within the cut-off procedure. The requirement of zero value of this 
remainder is then natural. It should be emphasized that both of these motifs are in 
accord with the statements discussed in the literature [HIE]- In other words, our latter 
treatment, albeit original, is developed in a quite traditional manner, without recourse 
to any additional adjustable function simulating the contribution of the lattice sum 
over direct space [8]. 

It is worth noting that according to (|65|) . the value of the spreading parameter 



is connected with the value of the cut-off parameter m. It is evident that this 

relation could be formally inverted, i.e., m might be treated as a function of a|™J as 
well. However, contrary to ar n -\ that is continuous initially, m is a discrete parameter 
and so a discrete set of values appropriate to different m actually arises as a 



solution of (|65|) . 

One more possibility to optimize the calculation of lattice series is associated 
with the fact that the lattice sum in direct space converges faster then that in 
reciprocal space. Therefore it may be advantageous to choose the cut-off parameter 
m R truncating the direct lattice sum larger than the cut-off parameter m applied to 
the corresponding reciprocal lattice sum. In this case the correlation between the 
cut-off parameter m R and the spreading parameter a/ n \ may be of interest, providing 
that the third parameter m is fixed [H [HI E] ■ The foregoing treatment based on (ffi5)) 
can be easily extended to the present case if we adopt 

r.[m+m a , m R +m a ] p[m, ma] _ n ( RR\ 

C (n) C (n) ~ U ' y Q0 > 

where we conjecture that the increment m a in both m and m R is the same, for 
simplicity. Note that the two limitimg cases mentioned above as associated with 
m a still take place. As far as the choice of m R is concerned, it turns out that this 
value is to be as large as possible. Actually its value is restricted by machine accuracy. 
On the other hand, its effect on ar n -\ is strong, as will be shown later on. 

In the events of confined spreading functions restricted by Rq , the Bertaut version 
dealing with non-overlapping charge distributions is the most popular. As is well 
known, the sum over direct space is then zero and so it does not contribute to the 
result that is determined by the sum over reciprocal space alone [T2l [Ti] . Keeping 
in mind that the value Rq must still be as large as possible, one may conclude that 
Rq is to be half the nearest interatomic distance pT2l [13l [16l [18l [20l [21l [28l [29l [30] . 
This convention is sustained in modern papers devoted to this subject as well [24]. 
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As a result, the assessment of the accuracy of computation is reduced to the classical 
consideration of the cut-off effect upon summation over reciprocal space [151 [T9"] . 

However, it is expedient to note that the aforementioned restriction on R may 
be convenient, but is not principal. Indeed, some further growth of Rq results in the 
appearance of a certain direct sum that is finite and can be calculated rigorously. 
On the other hand, the larger Rq the faster convergence of the series over reciprocal 
space. This is the reason to make use of the largest value of Rq consistent with machine 
accuracy of evaluation of the direct sum [2j [TBI El] • The latter is the problem typical 
of practical summation. 



5. Effect of infinite normalized spreading functions 



Here we apply the foregoing results to the classical NaCl structure composed of point 
charges, which remains attractive as a model system [21] 130] [31]. As known, a face- 
centred cubic Bravais lattice describes this structure, with two sites per unit cell. If 
d is the edge of an elementary cube, then v = rf 3 /4 and the basis vectors for point 
charges ±q may be chosen as: 

+ q: bi = (0,0,0), -q: 6 2 = (rf/2, 0, 0). (67) 
In terms of the elementary translations 



ai 



d(l,l,0) 



a 2 



d(l,0,l) 



a 3 



d(0,l,l) 



(68) 



2 ' 2 ' 

an arbitrary translation vector is of the form 

Ri = miai + m 2 a 2 + m 3 a 3 , (69) 

where rtij are integers. The elementary reciprocal lattice translations appropriate to 
the vectors in (|68|) are defined by the scalar product (o,ihj) = Sij, where Sij is the 
Kronecker delta, and are as follows: 

fc-^- <™> 

They compose a general reciprocal lattice vector of the form 

h = niihi + m 2 h2 + m^hz- (71) 
Note that definitions ([6"5]) and ([?0"|) are conventional [52] . 

On taking formulae ©, ^ and (p3T[) into account, expression (fTI)]) for the energy 
of interest is transformed into 



WW<\Ri\) W^(\Ri + b 2 \) 



\Ri 



\Ri 



q 2 ^ {l-cos[27T(hb 2 )}}S n (h) 2 

£ ^ = ™^ w q 

-q 2 & n \0). (72) 
According to (f55|) — (fTTj) . the parameters necessary for the further summation take 
the form 



\Ri\ = 

\R,+b 2 \ = 



(mi + mi) 2 + (m 2 + m 3 ) 2 + (m 3 + mi) 
d 



1/2 



1/2 



(mi + m 2 + l) 2 + (m 2 + m 3 ) 2 + (m 3 + mi) 2 
(mi + m 2 - m 3 ) 2 + (m 2 + m 3 — ?7ii) 2 + (m 3 + mi — m 2 ) 2 

1 - cos[2ir(hb 2 )} = { 



W-3 



1/2 



at mi + m 2 — m 3 even, 
at mi + m 2 — m 3 odd. 



(73) 
(74) 
(75) 

(76) 
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Table 1. The Coulomb energy £, in units of q 2 /d, per unit cell of the NaCl 
point-charge lattice. The calculation is based on a Gaussian charge spreading 
with = proposed originally by Ewald (£e w ) an d with fj, = •/ttv~ 1 ^ in 

agreement with our proposal (£our). The results, with significant figures only, are 
shown in dependence on m restricting actual ranges of summation over direct and 
reciprocal space, in accord with condition (1771 . 



m 




£*our 


1 


-3.5 


-3.5 


2 


-3.49513 


-3.4951292 


3 


-3.49512918927 


-3.49512918926636 a 



a This value agrees with the result of Sakamoto [33] . 



Note that the latter relation is very specific [21]. Employing relations (f73 |) -([76 |) 
in formula (|72p , we propose that the parameters of summation are restricted by a 
common condition 

\rrij | < to, (77) 

where an integer to is varying. The rate of convergence of series in equation (|72[) will 
be then studied in dependence on to. 

Here we start from a Gaussian spreading function that is a classical example 
[21 [13] . Bearing in mind that the effect of multiple charge spreading is reduced only 
to some definite scaling of Gaussian parameters in relations ([27]) -([29 ]) . we consider the 
case of a unique optimal value of a Gaussian parameter described by (|62[) . The value 
of the Madelung energy is shown in table [1] in dependence on m. The effect of the 
original value of /i = proposed by Ewald [H [7] is demonstrated in the same 

table for comparison. Table 1 shows that the rate of convergence is actually fantastic 
in both these cases, though our choice appears to be somewhat more advanced. Such 
results are in general expected [19j . This is the reason that we will not discuss further 
possibilities of optimization addressed to a Gaussian spreading function [S] . 

At this stage it is worth noting that the accuracy of results specified by Gaussian 
spreading function drops drastically if we restrict ourselves to spherical domains of 
summation over both the reciprocal and direct lattices. This fact is known [5]. 
Unfortunately, such a loss in accuracy caused by spherical modes of summation in the 
series at hand turns out to be typical of all the other cases discussed below. It means 
that modes of summation supporting the crystal symmetry are more advantageous 
for perfect crystals. On the other hand, very popular schemes of spherical summation 
seem to be essential in applications describing disordered systems Qj] . 

In the cases of a simple exponential spreading formulae ([20]) ([24]) and (|64p are 
utilized in ({72"]) first of all. The results of summation are listed in table [5] As 
anticipated, the tendency towards enhancing the rate of convergence takes place upon 
increasing the order of multiplicity n. This effect is quite general and has been pointed 
out by Bertaut [2] in connection with the particular cases of n — 1 and n — 2 at 
spatially confined spreading functions, as will be discussed in the sequel. 

However, table [5] also shows that, contrary to a Gaussian spreading function, the 
application of a simple exponential spreading with fixed optimal values of cx.i n ) results 
in much more moderate rate of convergence of lattice series in question. Therefore 
further efforts towards optimizing lattice calculations seem to be necessary here. The 
next step of optimization discussed is associated with a possible variation of «(„) in 
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Table 2. The Coulomb energy £. in units of q 2 /d, per unit cell of the NaCl point- 
charge lattice. The calculation is based on a simple exponential charge spreading 
with optimal values of ct( n \ described by relation (1<33 I t for n = 1, n = 2 and n = 3. 
The results are shown in dependence on m in condition l|77| l again. 



m 


£(i) 


£(2) 


£ (3) 


5 


-3.5 


-3.4951 


-3.495129 


10 


-3.5 


-3.495129 


-3.495129189 


15 


-3.495 


-3.4951292 


-3.49512918927 


20 


-3.495 


-3.49512919 


-3.4951291892664 


25 


-3.495 


-3.495129189 


-3.4951291892664 


30 


-3.4951 


-3.4951291893 


-3.4951291892664 


35 


-3.4951 


-3.4951291893 


-3.49512918926636 a 



a This value agrees with the result of Sakamoto |33| . 



dependence on m. With the help of condition (f65|) . this effect can be readily taken into 
account. Here two particular cases of m a = 1 and m a = 10 are studied. The results 
of calculation are compiled in tabled We see that the tendency towards reducing the 
value of oe/ n ) upon growing m is common to all the events of n under consideration. 
It implies that the effective charge distributions become more diffuse as the range 
of summation increases. On the other hand, the multiple spreading leads to larger 
values of ctr n \. This output agrees with the results of It is important that the 

enhancement of the rate of convergence arises in both the cases of m a , but the effect 
at m a = 10 is somewhat stronger than that at m a = 1. 

The latter trend can be fruitful when we go over to the next step of optimization 
based on relation ([66]) . The event of m — 1 and m R = 22 is studied, providing that 



Table 3. The Coulomb energy £, in units of q 2 /d, per unit cell of the NaCl 
point-charge lattice in the case similar to that in table [2] but with c*( n ) of which 
variation with m is specified by (165 I t at m a = 1 and m a = 10, respectively. 



m 




£ W 


«(2) 


£ (2) 


"(3) 


£ (3) 


case of m a = 1 


1 


6.8265 


-3.5 


11.2028 


-3.5 


14.4677 


-3.495 


3 


5.0228 


-3.5 


8.7002 


-3.49513 


11.5983 


-3.4951292 


5 


4.0395 


-3.495 


7.2189 


-3.4951292 


9.8319 


-3.49512919 


7 


3.4073 


-3.495 


6.2024 


-3.49512919 


8.5572 


-3.49512918927 


9 


2.9639 


-3.49513 


5.4620 


-3.49512919 


7.6008 


-3.495129189266 


11 


2.6336 


-3.49513 


4.8970 


-3.4951291893 


6.8568 


-3.4951291892664 


13 


2.3769 


-3.49513 


4.4502 


-3.49512918927 


6.2618 


-3.49512918926636 a 


case of m a = 10 


1 


6.5898 


-3.495 


11.1262 


-3.4951292 


14.4308 


-3.4951291893 


2 


5.5308 


-3.495 


9.6293 


-3.4951292 


12.6762 


-3.49512918927 


3 


4.8102 


-3.495 


8.6223 


-3.49512919 


11.5650 


-3.49512918927 


4 


4.2661 


-3.4951 


7.7974 


-3.49512919 


10.5979 


-3.49512918927 


5 


3.8453 


-3.4951 


7.1272 


-3.495129189 


9.7854 


-3.4951291892664 


6 


3.5093 


-3.49513 


6.5721 


-3.4951291893 


9.0943 


-3.4951291892664 


7 


3.2346 


-3.49513 


6.1058 


-3.4951291893 


8.5019 


-3.4951291892664 


8 


3.0052 


-3.49513 


5.7084 


-3.4951291893 


7.9888 


-3.49512918926636 a 



a This value agrees with the result of Sakamoto ,33; ■ 
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Table 4. The Coulomb energy £, in units of q 2 /d, per unit cell of the NaCl 
point-charge lattice in the case of a simple exponential charge spreading like that 
in table [3] but with m = 1 and m R = 22 restricting lattice summations over 
reciprocal and direct space, respectively, in accord with (|77j l. Then the values of 
cti n \ at n = 1, 2 and 3 are determined by II66I I at m a = 1 and m a = 10. 





n 






1 


1 


1.0346 


-3.4951 




2 


1.9312 


-3.4951292 




3 


2.6808 


-3.495129189 


10 


1 


1.0681 


-3.49512919 




2 


1.9562 


-3.49512918927 




3 


2.6986 


-3.495129189266364 a 



a This value agrees with the result of Sakamoto 33 . 



the values m a = 1 and m a = 10 are utilized in ()66j) . The corresponding results are 
listed in tabled] Table [4] shows that the values of a(„) corresponding to each value of 
n appear to be close, but the effect of m a = 10 is more pronounced again. 

6. Trends in the application of confined spreading functions 

As far as confined spreading functions are concerned, here we restrict ourselves to 
the situation opposite to the classical one described by charge distributions non- 
overlapping after spreading [121 1 1 4] . In other words, we propose that the parameter 
i?o in (|30p is greater than the lattice parameter d. It is important that the sum 
over direct space is always finite in such a case and so it can be counted with a high 
precision. In practice our choice of Rq is actually restricted by machine accuracy. As a 
result, optimum values of Rq appear to be different and depend on each particular case 
under consideration. For example, if k = 0, then we adopt Rq = Id upon considering 
all three cases for n from one to three. Note that the value of m in (fTT)) restricts 
only the reciprocal lattice series now. The computation based on formulae (|3Tj) - (|40)) 
substituted into (|72[) gives rise to the results presented in table Table shows that 



Table 5. The specific Coulomb energy £, in units of q 2 /d, for the NaCl point- 
charge lattice is obtained in dependence on the restricting parameter m at a fixed 
value of i?o = 7d common to all the cases of the confined polynomial spreading 
function go(2)( r ) defined by formula l|31[ l. The cases of £( n ) at n = 1, n = 2 and 
n = 3, which are specified by (132 1 — (140 l l in formula (1721 are considered. 



m 




%) 


£ (3) 





-3.5 


-3.49513 


-3.495129189 


1 


-3.495 


-3.4951292 


-3.49512918927 


2 


-3.495 


-3.4951292 


-3.49512918927 


3 


-3.495 


-3.4951292 


-3.4951291892664 


1 


-3.495 


-3.49512919 


-3.4951291892664 


5 


-3.495 


-3.49512919 


-3.4951291892664 


10 


-3.4951 


-3.49512919 


-3.495129189266364 


15 


-3.495129 


-3.4951291893 


-3.4951291892663644 a 



a This value agrees with the result of Sakamoto [33_ (see also references |34l 1351 ). 
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Table 6. The specific Coulomb energy £, in units of q 2 /d, for the NaCl point- 
charge lattice obtained with making use of either 91(3) (r) defined by H41H at 
Ro = 6.5d or 92(3) ( r ) denned by (149 I t at Ho = 7d, in dependence on the cut-off 
parameter m. The cases of n = 1 and n = 2 for both these spreading functions 
are considered together for comparison. 



in 


case 


of Sl(3) 0) 


case 


of 92(3) 0) 




£(2) 


%) 


£(2) 





-3.495 


-3.4951292 


-3.495 


-3.4951292 


1 


-3.49513 


-3.4951291893 


-3.4951 


-3.495129189 


2 


-3.49513 


-3.4951291893 


-3.49513 


-3.49512918927 


3 


-3.495129 


-3.49512918927 


-3.49513 


-3.49512918927 


1 


-3.495129 


-3.49512918927 


-3.49513 


-3.49512918927 


5 


-3.495129 


-3.49512918927 


-3.49513 


-3.49512918927 


10 


-3.495129 


-3.4951291892664 


-3.495129 


-3.4951291892664 


15 


-3.4951292 


-3.4951291892664 


-3.495129 


-3.4951291892664 


20 


-3.49512919 


-3.49512918926636 a 


-3.49512919 


-3.49512918926636 a 



a This value agrees with the result of Sakamoto 33 



the accuracy of the output obtained in these cases enhances rapidly upon changing n 
from 1 to 3. This result agrees with the general tendency mentioned above [Mj. By 
the way, it is curious that the value at the bottom of the first column appears to be 
in between the values determined by to = and to = 1 in the second column and the 
same is right for the second and third columns as well. Of course, this fact seems to 
be occasional, but it is expressive. 

Note that the case of m = here and below describes the contribution of the sum 
over direct space alone. The corresponding results shown in table [5] give evidence that 
the summation over real space alone can be very precise and so may be treated as one 
more efficient approach to the direct lattice summation problem [H [T31 [35] . Within 
such a treatment the series over reciprocal lattice vectors may in turn be regarded 
as an additional correcting contribution to the direct-space one. It is natural that it 
leads to the further essential enhancement of the overall accuracy. 

At the next step we consider the same energy counted with making use of the 
spreading function 31(3) (r) specified by ([41]) and (|42|) at R = 6.5d as an optimal value. 
The corresponding cases of n = 1 and n = 2 are determined by formulae ([4"51) - (|4"4")) 
and ([4"5|) - (|48]) . respectively. We also examine the case of 32(3) ( r ) that is optimized at 
Ro = Id and is determined by relations (f49|) - ([56|) . with including the events of n = 1 
and n = 2 as well. The results of calculation are given in table El On comparing these 
results with the first two columns of table it is evident that the latter issues are 
more efficient in the case of n — 1 and the same is right for n = 2. The tendency that 
every case of n = 2 is much more accurate than the corresponding case of n = 1 is 
also maintained. On the other hand, table[6]shows that any advantage of 32(3) ( r ) over 
ffi(3)( r )i that is expected for non-overlapping spreading functions [24], is not realized in 
the cases of large Rq. Moreover, the case of 31(3) (r) at n — 2 appears to be somewhat 
more precise, keeping in mind that this effect is achieved at a less value of Rq. 

It is significant that the enhancement of accuracy takes place along the set of 
spreading functions considered in |Appcndix A| even at n = 1. Note that the formulae 
associated especially with (74(5) (r) and 3.5(6) (r) are rather novel. This is the reason 
that the results appropriate to those cases and obtained at the corresponding optimum 
value of R = 15d are listed in table [7] Indeed, the rate of convergence increases as 
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Table 7. The specific Coulomb energy £, in units of q 2 /d, for the NaCl point- 
charge lattice is obtained in dependence on the restricting parameter m at a fixed 
value of Rq = 15d common to all the cases of the confined polynomials specified 



by 


I|A.1||-||A.12|| in|Appcndix 


A|at n = 1. 




m 










case of g 2{4) (r) 


case of g 4{5) (r) 


case of S 5(6 )(r) 





-3.49513 


-3.49512919 


-3.49512919 


1 


-3.495129 


-3.49512919 


-3.495129189 


2 


-3.4951292 


-3.4951291893 


-3.4951291893 


3 


-3.4951292 


-3.4951291893 


-3.49512918927 


4 


-3.49512919 


-3.49512918927 


-3.49512918927 


5 


-3.49512919 


-3.49512918927 


-3.49512918927 


10 


-3.49512919 


-3.4951291892664 


-3.4951291892664 


15 


-3.495129189 


-3.49512918926636 


-3.4951291892664 


20 


-3.495129189 


-3.49512918926636 


-3.49512918926636 


25 


-3.495129189 


-3.49512918926636 


-3.49512918926636 


30 


-3.495129189 


-3.49512918926636 


-3.495129189266364 


35 


-3.4951291893 


-3.49512918926636 


-3.4951291892663644 a 



a This value agrees with the result of Sakamoto |33| (see also references 1341 1351 ). 



k grows. Of course, this effect turns out to be less prominent in comparison with 
the cases of n = 2 and n = 3 mentioned above. Nevertheless, at k — 5 the limiting 
precision adopted in our calculations is eventually attained at as well. 

Note that if n is odd, then the contribution of the reciprocal lattice sum tends 
to the exact value in an oscillatory manner due to trigonometric functions describing 
Sfc(s) (A) appropriate to all the events in this section [T3J HI] ■ 

Some final comment is necessary on the question why all the foregoing tables 
contain limiting energy values with a rather large number of significant figures. 
Indeed, one may think that the lattice parameters of real structures are usually 
known up to four, at best six, figures only [37]. Nevertheless, even in this case the 
numerical calculation should be a bit more accurate. However, the reason of our 
ultimate accuracy is quite different. As shown, the analytic accuracy of a series is 
associated with the number of unit cells taken into account and restricted by a cut- 
off parameter introduced upon series computation. On the other hand, restricted by 
machine accuracy, the overall accuracy of computation depends on the total number 
of implemented operations and so is connected anyhow with the total number of point 
charges taken into account. The interplay between these tendencies implies that the 
effect of machine accuracy may be predominant if the number of charges per unit cell 
increases. In other words, an ultimate accuracy achieved for model structure with a 
small unit cell may be regarded as a guaranty of a sufficient accuracy while each of the 
approaches developed above is applied to modern compounds with large unit cells. 

7. Conclusion 

In summary, it has been confirmed that the effect of a multiple charge spreading results 
in increasing the rate of convergence of the lattice series. In other words, we deal with 
one more way to make the convergence faster, providing that this effect becomes 
stronger and stronger with every further repetition of the procedure of spreading. It 
is important that this is a property common to both the classes of infinite and confined 
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normalized spreading functions. 

The optimization of the shape of spreading is more traditional in the problem 
of enhancing the rate of convergence. Nevertheless, we have proposed some novel 
approaches to this task. In particular, for infinite normalized functions of spreading 
two different situations are considered. In the case of a fixed spreading parameter, the 
universal approach of optimization is based on the investigation of the concomitant 
terms contributing to the Coulomb energy, which are independent of the lattice sums. 
The problem of optimizing the spreading parameter in dependence on the cut-off 
parameters of summation has been solved with the help of some conditions imposed 
on the remainder of the Coulomb energy, providing that this remainder is also treated 
in a truncated form. In the case of confined spreading functions the optimization in 
question is reduced to separating any main set of polynomials ensuring the progressive 
enhancement of the rate of convergence of the sum over reciprocal space. 

It is found that the effect of optimization becomes much more efficient if the 
multiple spreading is applied as well. The only case independent of the multiple 
spreading is described by a Gaussian spreading function due to its invariance with 
respect to spreading in a multiple manner .22] . Moreover, it turns out that the 
convergence with a Gaussian spreading function is the most prominent even at a fixed 
optimized spreading parameter. This fact gives evidence that a Gaussian spreading 
seems to be the most suitable one in the treatment based on the charge spreading as 
a whole. 

On the other hand, the investigation of confined functions of spreading shows that 
it is fruitful to choose larger values of the parameter restricting those functions. Of 
course, a certain finite sum over sites of the direct lattice then appears. However, there 
is no problem with its calculation. On the other hand, the contribution of the sum 
over reciprocal lattice becomes smaller, up to the case where the latter sum describes 
only a small correction to the result of direct summation. Note that this trend is 
eventually common to all the examples discussed above. By the way, it implies one 
more approach to the problem of direct summation of Coulomb series in crystals. 

Appendix A. Particular forms of a polynomial spreading function 

Here we continue the fundamental series of polynomials beginning with (|31[) and (141[) , 
which produce S(h) with the enhancing power of Y in the denominator, providing 
that the polynomial degree be a minimum. Note that the case of n — 1 alone will 
be considered here. So, the next example of such a polynomial, with including the 
concomitant relations, is of the form 



52(4) (r) 




15(1 -x) 2 

2tzRI 
60 / 



3 



Y 



(A.l) 



(A.2) 




at < z < 1, 



at z > 1, 



(A.3) 



0> = db 



(A.4) 



where the notations of section [3] are utilized. 
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The next polynomial of spreading within the set under consideration is not simply 
defined by (1 — x) 3 , but it is described as 

105(1 -x) 3 (l + 3x) 



94(5) (r) 



WirR 3 



with the concomitant quantities 

630 /_. 
Y5\ 



^4(5) (h) 



( sin Y - 
(l-zf( 



8 + 7cosy 15 sin Y 



Y Y 2 
19z + 15z 2 ) 





ad) (o) - 2L 



at < z < 1, 
at z > 1, 



The relations specifying one more pattern of this set are of the form 

21(1 -x) 4 (l + 4.x) 



05(6)0) = 



'5(6) 



(h) 



5040 r 



ye 



2nR 3 
4 - cos Y 



9sinF 24(1- cos Y) 



Y 



5(6) 



(l-z) 6 (l + 3z + 32 2 ) 




Y 2 

at < z < 1, 
at z > 1, 



5(6) V / ^ 



(A.5) 

(A.6) 
(A.7) 
(A.8) 

(A.9) 
(A.10) 
(ATI) 
(A.12) 



Interested in the cases of s < 6, we will not propose the complete set of values 
characteristic of the next example. However, it is worth pointing out that the structure 
of the corresponding spreading function becomes of a more complicated form again 
and is as follows: 

45(l-x) 5 (l + 5x + 8a; 2 ) 



57(7)0) 



471-i? 3 , 



(A.13) 



The corresponding concomitant quantities take the form 

75600 



SV( 7 ) (h) 



Y 7 

192(1- cos Y) 



24 -15 cos Y 87 sin Y 
sm Y H — h 



Y 



Y 2 



Y 3 



[jW (0) - J^t 



(A.14) 
(A.15) 



and are sufficient to describe point-charge lattices at least while different spreading 
functions of form (|A.13[) do not overlap. 
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